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Current methods to describe the thermodynamic behavior of many-particle systems are often 
based on perturbation theory with an unperturbed system consisting of free particles. Therefore, 
only a few methods are able to describe both strongly and weakly correlated systems along the same 
lines. In this article we propose a cumulant approach which allows for the evaluation of excitation 
energies and is especially appropriate to account for the thermodynamics at low temperatures. The 
method is an extension of a cumulant formalism which was recently proposed to study statical and 
dynamical properties of many-body systems at zero temperature. The present approach merges into 
the former one for vanishing temperature. As an application we investigate the thermodynamics of 
the hole-doped antiferromagnetic phase in high-temperature superconductors in the framework of 
the anisotropic t-J model. 
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I. INTRODUCTION 

For the investigation of many-particle systems at finite temperatures the free energy F plays a central role. One 
basic property of the free energy is its size consistency, i.e. the free energy scales with the size of the system. Each 
approximation which is used to evaluate F must preserve this property. In diagrammatic approaches size consistency 
is guaranteed by the fact that in any physical quantity only linked diagrams enter. Usually, diagram technique makes 
use of Wick's theorem, which is appropriate only if the unperturbed Hamiltonian is a single-particle Hamiltonian. 
Therefore standard diagrammatic approaches are restricted to weakly correlated systems. -Aa alternative approach 
to evaluate statical and dynamical properties at zero temperature was recently proposedou and is based on the 
introduction of cumulants. As is well known from statistical physics the use of cumulants ensures size consistency. 
This favors the cumulant method especially for the description of strongly correlated systems though it may be applied 
as well to weakly correlated systems. The main aspect of the present paper is to investigate excitation energies which 
are needed to evaluate the free energy. 

There is a natural danger that a cumulant expansion of F is only valid for high temperatures. Expanding with 
respect to small values of the inverse temperature (3 = 1/kT the partition function reads 

Z = Z(0)(l-l3H+^ 2 H 2 ---)o , (1) 

where (• • -) = Tr(. . -)/Tr 1 stands for an unweighted average at infinite temperature and Z (0) = Tr 1 is the trace 
of the unity operator 1. Expanding \nZ, suitably regrouped into contributions in powers of /?, one finds for the free 
energy 

F = ~ In Z (0) + (H) + |s [(H\ (2) 
+ ^ [(H 3 ) - 3(H 2 ) (H)o + 2(H 3 )o] +■■■ , 

Note that the expressions in the brackets [• • •] are cumulants. 

In the present approach we want to avoid that a free energy expansion like (0) is only valid for high temperatures. 
Of more interest are low-temperature expansions. A low-temperature expansion for F must reduce to the ground- 
state energy in the limit T — > 0. In contrast to a high-temperature expansion where (3 is small, for sufficiently low 
temperatures the canonical weight exp (—(3(E n — Eq)) is small for excitations E n higher than the ground-state energy 
Eq. Therefore, in case of a non-degenerate ground state the free energy can be expanded for low temperatures as 
follows 

F = -- In Z = E - - In (1 + £ e -/><s»-*0) (3) 

P P n>0 
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= ^-4(E e "" (M) +-) 

P n>0 

A neglection of high-lying excitations is not always allowed. For instance, for the one-dimensional Ising model in a 
longitudinal field the thermodynamic limit and the limit T — > do not commute. This is an alternative description 
of the fact that long-range order is destroyed in one dimension. 

In this paper we propose a cumulant formalism for the calculation of excitation energies of correlated electronic 
systems. This method is based on a perturbational approach, i.e., the Hamiltonian is splitted into H and H\ 
with H being exactly solvable. One starts from eigenstates of H and includes the effect of Hi by an exponential 
ansatz. Especially for the zero-temperature version of the method this can be shown to be equivalent to summing 
a perturbation series to infinite order. So the present method is well suited for systems which can in principle be 
treated perturbatively, but it can not account for systems with very large fluctuations, e.g., most one-dimensional 
spin models. 

This paper is organized as follows. In the Sec. II we shortly review the size-consistent ground-state version of the 
cumulant method. In Sec. Ill the cumulant method for the computation of excitation energies is presented. With 
this general scheme we are able to calculate partition functions. To demonstrate the applicability of the method 
we investigate the t-J model at weak doping as a topic of current interest. We consider a two-dimensional model 
with anisotropic magnetic exchange and calculate the staggered magnetization within the antiferromagnetic phase in 
dependence on temperature and hole concentration. A discussion and concluding remarks are put in the last section. 



II. REVIEW OF THE ZERO-TEMPERATURE VERSION OF THE CUMULANT METHOD 



a better understanding of the cumulant method its ground-state version is briefly reviewed. For more details 
. The method starts from the definition of the function 

/(A)=ln(0o|e- A(ffo+ffl) e AHo |0 o ) (4) 

where Hq is the rigorously solvable unperturbed Hamiltonian Hq with the ground state \4>o), i-e., i?o|0o) = £o|0o)- 
The aim is to calculate the ground-state energy Eq of the full system, i?|V>o) = I V'o) ■ The shift of the ground-state 
energy 5Eq = Eq — eo due to the perturbation Hi can be derived in a straightforward way. Introducing the Liouvillian 
Lq which is defined by LqA — [Lq, A] for any operator A, equation (^) is transformed into 

/(A)=ln(0 o |e- A(Hl+io) |0o> . (5) 
Next, we define the Laplace transform of the function / (A) by: 

/>) = - / /(A) e Xz dz , K{z}<0 . (6) 



o 



One can showEJu that the energy shift 5 Eg with respect to the unperturbed ground-state energy eo is given by 

5E = lira z 2 f » . (7) 

z— >0 

On the other side, equation (|^) is used to express 5Eq in terms of cumulants 

SE = lim ^ \H! (l + -i -Hi) |0 O ) C . (8) 



2^0 \ z — Hi — Lq 

Here (</>o|...|<^o) c denotes cumulant expection values with respect to the unperturbed ground state \4>q). Cumulant 
expectation valuesB for a product of arbitrary operators Ai with an arbitrary state are defined by: 

(^XlArW = fn(^) n jln^ine^^)U i= o Vi . (9) 

For a note on generalized cumulants see Appendix A. The quantity inside the bracket of (||) is called wave operator 
fi, (it has similarity with the Moller operator known from scattering theory), 



2 



= 1 + lim 



2^0 z — Hi — L 



(10) 



Thus we can rewrite SEq as 



SE = (&,|iJifi|<fo 



or E = {<po \Hn\(f> 



(11) 



Within cumulants, the operator f2 transforms the ground state \(/)q) of the unperturbed Hamiltonian H into the 
full ground state \ipo) of H . Expanding ( ^p| ) into powers of H\ it can-,be shown that ( JTT| ) is equivalent to Rayleigh- 
Schrodinger perturbation theory summed up to infinite order, see e.g.u. 

There is no general rule how to split H into Hq and Hi except that the overlap between the unperturbed and the full 
ground state has to be non-zero, i.e., (ipo\(f>o) ^ 0. The operator £1 describes the influence of Hi onto \<fio). This effect 
should be a small correction to |0o), i-e., it should be treatable perturbatively in the sense that it can be obtained 
by summation of a perturbation expansion. This is usually fulfilled if Hq is the dominant part of the Hamiltonian. 
Therefore, for strongly correlated systems Ho should consist of the correlations (or at least part of them) whereas Hi 
usually contains the hybridization. In the application presented in Sec. IV dealing with the weakly doped t — J model 
Hq contains the Ising part of the magnetic interaction whereas Hi consists of its transverse part and the electron 
hopping. Other interesting examples may be found in refs.Qu. If the full ground state \ipo) is expected to break a 
symmetry of H then Ho has to be chosen so that its ground state |0o) breaks this symmetry, too. n 

Instead of using the explicit form ( |l0| ) of the wave operator Q an exponential ansatz was proposedu 



n 



s 



(12) 



where {S^} is a set of relevant operators. They have to be chosen in such a way that exp(J^ 



a 



(with 



appropriate parameters a^) represents a good approximation of the exact ground state. The yet unknown parameters 
an, are to be determined from the following set of equations 



bo\StHn\. 



= 0, i/= 1,2,3, 



(13) 



These equations follow from the condition of £l\(f>o) being an eigenstate of H. Note that equations (Jll|) and ( |l3| ) allow 
for the computation of the ground-state energy. For algebraic reasons it is suitable to use operators which only 
create (but do not destroy) fluctuations with respect to the unperturbed ground state \4>o). Then the expansion of the 
exponential in ( |ll| ) and ( |l3[ ) stops after a few terms because the Hamiltonian H and the adjoint fluctuation operators 
Si have to remove all fluctuations which are created by O and Hi in the state |</>o). 

The choice of appropriate operators S v is most important for actual calculations using the cumulant method. These 
operators describe fluctuations introduced into \4>o) by successive application of Hi. In principle, they can be derived 
systematically from the explicit form ([l0]) of the wave operator Q. For practical applications this might be only of 
little help especially if the main physical effect comes from higher powers of Hi. In such a case a small set of few 
relevant operators leads to a far simpler description of the main effect than including a large set of powers if™. The 
selection of relevant operators for the cumulant method can be seen similar to the selection of dynamical variables for 
Mori-Zwanzig projection technique: Formally, variables can be systematically derived from the Liouvillian, but often 
choosing variables from physical insight is more useful. 

For a discussion of the cumulant method and its relation to other methods like variational calculations and the 
coupled-cluster method see appendis^ the past, equations (|ll| , |l3|) was used to evaluate ground-state properties 
of several systems, see for instanceaE 



III. EXTENSION TO FINITE TEMPERATURES 



Now we present a cumulant scheme for calculating excitation energies. To discuss the influence of the perturbation 
Hi we formally introduce an additional parameter A in the Hamiltonian: 

H = H + XHi (14) 

As is shown below, the wave operator f2 from the zero-temperature approach will be replaced by a unitary operator 
U which diagonalizes the full Hamiltonian H. The operator U transforms all eigenstates of the unperturbed system 
into eigenstates of the full Hamiltonian. Being |0„) and \ip n ) the eigenstates of Hq and H, 

H Q \(j) n ) =e n \(j) n ) , H\r/j n ) = E n \r/j n ) , (15) 
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the action of U is defined by 



- u\<t> n ) 

with U depending smoothly on A and U — > 1 for A — > 0. Every unitary operator can be written as 

U = e s with 5 f = -5. 



(16) 



(17) 



In general, the unitary transformation U is not known for a given system. Therefore approximations for U have to 
be used. Generalizing the exponential form (fL2|) for the wave operator f2 from the zero-temperature approach to a 
unitary operator, we make the following ansatz 



U 



(18) 



where both 5^ and the adjoint operators Si are enclosed in the exponential. The {S v } form a set of relevant operators 
as in the zero-temperature approach. In most cases, using a finite number of operators S v represents an approximation 
of U. However, with the number of operators going to infinity the exact transformation U can be approximated with 
arbitrary accuracy. 

In contrast to the ground-state approach the operators and can create and also annihilate fluctuations. This 
also means that the application of (S^ — SJJ to an eigenstate \4> n ) of Hq (with energy e n ) leads to a new state which 
may overlap with eigenstates of Ho having energies both below and above e„. 

In the finite-temperature approach the former equations (JTA) and ([L3j) are replaced by 



E n = (HU) c n , 

o = ((s„ - stymrx 



v= 1,2,3,. 



(19) 
(20) 



where (• • -) n refers to the cumulant average now formed with the unperturbed eigenstate \4> n ) of Ho, i.e., (• ■ -) n = 
(4> n \ ■ ■ ■ \4>n) c - The parameters a n can be determined from Eqs. (|2(i|). Note that the cumulants in ( pfj| ) can be formed 
with any unperturbed eigenstate \<p n ) of Ho as long as the exact expression for the unitary transformation U is known. 
However, for any approximation for U the parameters may vary with a different choice of the \4>n) used for the 
cumulant expectation values. I— ■ 

To prove Eqs. (JT^j2FJ) one can either make use of a method recently proposedlij which is based on integrating over 
infinitesimal transformations (1 + S/N). Here we rather prefer to make explicit use of the definition (0) of cumulant 



expectation values. Expanding U 
one finds for the r.h.s. of ( |l9|) 



in powers of S and taking the definition <M) for cumulant expectation values 



^ 711 I 



(21) 



1 



lim 



d d 



in 



ln(0„|e AlH e x * s \<f> n 



Summing over m one immediately finds the desired result 

d 



{He 



S\c 



d 

d 
9AT 



In > 
In > 
In. 



p XiH e (A 2 + l)S| 



(22) 



A 1= A 2 =0 



Ai=0 

= E n 



A,=0 



The second equation (|20|) can be proved in close analogy. 

The changes imposed by the replacement of 51 by U give rise to some problems for actual calculations. The expansion 
of U does not terminate any longer after a few steps because the powers of SI in the unitary operator U can remove 
fluctuations created before by S^. For this reason an infinite series occurs. It may be possible, however, to find a 
solution in two steps: (i) First, the cumulants can be rewritten in terms of usual expection values formed with \4> n ) 
(see Appendix B), e.g. 
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n _ ,o rr ssc _ (SuHe s ) n (S„e s ) n (He s ) n 

\i-V v He) n - {eS)n {eS)n {eS)n [26) 

E n = {He s r n = , (24) 

\ e In 

with U = e s and S — ^ a^S^ — Sj). These relations can be easily verified by use of equation @. (ii) In a second 

step the exponential e s should be decomposed. The basic idea is to extract a factor next to \4> n ) that only annihilates 
fluctuations when it is applied to \<j) n ). This means that the power series of this factor usually stops after a few terms. 
As an example we consider the factorization for the special case of being a boson operator, [S^, S^] = 1. In this 
case a well known decomposition for U reads 

U = e Q (^- s i) = e ~ aS l e - ^ e aS " . (25) 

If U is applied to an eigenstate of Hq containing a finite number n of bosons, (SJJ 71 , the power series of the factor 
e M stops after n + 1 steps. For n = (ground state) only the identity operator would survive. Proceeding this way 
the application of U is much easier to evaluate than before. In general, decompositions like ( p5| ) depend on the details 
of the given problem. The example in Sec. IV will show that the outlined evaluation scheme is an appropriate tool 
to account for the influence of low- lying excitations especially for low temperatures. 

If we are faced with degenerate eigenstates of Hq with a degeneracy lifted by Hi we have the freedom to choose 
initial states \4> n ) for the cumulant method (|l9 20) from the subspace of degenerate states. In principle for every 



choice of \4> n ) a unitary transformation diagonalizing the Hamiltonian can be found. But if we use a certain ansatz for 
the operator U then the best choice of \c/> n ) depends on the form of this ansatz. Assume that the states \4> ni ), |<An m ) 
are degenerate eigenstates of H with energy e„. The correct states \<j) ni ) with U\(j) ni ) being eigenstates of H (for a 
certain form of U) are linear combinations of the \(f> ni )' 

m 

\K) = ^rfij\K)- ( 26 ) 

To find the-a'y one uses again Eq. ( po|) provided by the cumulant method. Note that Eq. ( |20| ) also holds for generalized 
cumulantstL3 (see Appendix A), i.e., defined with a bra vector different from the ket \4>n}' 

0= {<S>\{S V ~ HU\<p n ) c . (27) 

The arbitrary vector ($| should have a non-zero overlap with \<j> n ). Evaluating the cumulants in analogy to (23|24|) 
and using the linear combination ( p6| ) instead of the ket vector \4> n ) we obtain 

Y,Jv(®\S„HU\q> nj ) = E ni J2lij(®\S»U\4> nj ) ■ (28) 



For fixed U = e and appropriate operators £„, Eq. (|2£|) is a generalized eigenvalue problem for the E ni and 7^ and 
can be solved by standard methods. 

We should like to mention that the method presented here for calculating excitation energies on the basis of 
cumulants is rather different from the method proposed recently by Schork and FuldeuJ. There the zero-temperature 
wave operator Q ( jToj ) was applied to excited eigenstates \4> n ) of Hq to obtain the full eigenstate \i/) n ). This approach 
is only valid if \(il> n \(j) n )\ 2 > h which is a significant restriction compared to the present approach based on the 
introduction of the transformation U. 



IV. APPLICATION TO THE T-J MODEL AT WEAK DOPING 



The thermodynamics of the antifcrromagnetic phase in high-temperature superconductors is discussed as the major 
application of the present cumulant approach for finite temperatures. It is widely accepted that the electronic 
properties of these systems are mainly determined by the charge carriers in the CuC>2 planes—Of special interest is the 
interplay between antiferromagnetism and hole doping. Neutron scattering experiments! 1 3c3 show that the effective 
magnetization strongly depends on the hole concentration S within the Cu02 planes. Both the Neel temperature and 
the staggered magnetization decrease rapidly with increasing S and vanish at a critisaLconcentration S c . For reviews 
on experimental and theoretical investigations of high-T c superconductors see refs.tZHla. 
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The essential aspects of the low-energetic electronk.-d|3grees of freedom of the CuC>2 planes are by now believed to 
be well described by the two-dimensional t-J modeli3'E3: 

h = -t x; (*u-* + 4m + j e - ! r l ) c 29 ) 

As above, Si is the electronic spin operator and the electron number operator at site i. Note that the Hamiltonian 
(p9|) is defined in the subspace of the unitary space without double occupations of sites. The electronic creation 
operators c\ a are not usual fermion operators but rather exclude double occupancies: 

4 = 4(i-"i.-«o ( 3 °) 

At half filling the t-~L Hamiltonian reduces to the antiferromagnetic Heisenberg model. 

In a recent papertil we have investigated the behavior of the ground-state sublattice magnetization with increasing 
hole concentration using the cumulant formalism for zero temperature. Here we want to generalize these results and 
calculate the temperature and doping dependence of the magnetization. We consider—the situation of weak doping, 
i.e., of a few non-interacting holes present in the half- filled system. As done in ref.til we will use approximations 
equivalent to linear spin- wave theory for. expectation values with spinpperators and assume independent hole motion. 
Several analytical calculations at zeroE^tifLij and finite temperatures show that the strong coupling between holes 
and spin waves leads to a softening of long- wavelength spin excitations. This softening causes a reduction of the order 
parameter at zero temperature as well as a decrease of the Neel temperature with doping. 

Note that a non-zero spontaneous magnetizatiop,in a two-dimensional isotropic system at finite temperature is 
strictly forbidden by the Mermin- Wagner theoremca. The high-T c materials show a spontaneous magnetization at 
T > because of the anisotropy due to the magnetic coupling between the copper-oxide layers. For a simple description 
of that fact we introduce an anisotropy in the Heisenberg exchange, i.e., we use an exchange constant of J(l + e) 
in z direction. Within linear spin-wave theory this is equivalent to an additional staggered field eJ parallel to the 
z-axis. The value of the anisotropy e can _he_fitted from the Neel temperature Tjv observed in experiments. It can 
also be extracted from NMIL-data, see e.g.Ej'El Experimental values are J ~ Q.lbeV and Tjv ~ 400-ftT « 0.25J in the 
undoped materials (YBCO)ll3, so we are interested in the temperature range 0...0.5J. 



A. Undoped case 

First we briefly consider the half-filled case, i.e., the undoped system. The t-J model reduces to the s — 1/2 
Heisenberg antiferromagnet. For the application of the cumulant formalism we decompose the anisotropic Heisenberg 
Hamiltonian as follows: 

Hq = Hising 

= J(l + e) -1^), (31) 

<ij> 

Hi= H ± = J - J^i^S/ + S+Sr). 

<ij> 



We will use linear spin- wave approximation for the spin operators to calculate expectation values, i.e., the Fourier- 

S ± 



transformed spin operators are replaced by bosons a q and 6 q . These operators are spin- wave operators for the 
T-and j-sublattice defined by 



a+ = -jL= ye l ^S~ , (32) 

q tti 

b+ = — L= ws+ . 

q u 

In the approximation of linear spin-wave theory they obey usual Bose commutation relations: 

K><]- = <W [b^b+J- = <5 qiq2 , [a&>,6«]_ = 0. (33) 
Within linear spin-wave approximation the Hamiltonian ( |3l| ) takes the form 



G 



J (1 + e) z / A ^ , + + v \ 

#o = 2 ~ ^ + q q ^ J ' 

Hi = E7(q) M-q + <&t q ) (34) 
q 

where zo7(q) = 2 (cos + cos q y ) is the Fourier-transformed exchange coupling and zq the number of nearest neighbor 
sites. The wave vectors q have to be taken from the magnetic Brillouin zone. 

The ground state \4>o) of the unperturbed Hamiltonian H is the antiferromagnetically ordered Neel state. The 
excited states of Hq are obtained by creating fixed numbers of bosons a q and 6+ in the Neel state: 

l«W«*>> = II 4=T -7=7 ( 35 ) 

where {n aq ,n;, q } is the set of numbers of spin-wave excitations with momenta q on the |-and J,-sublattices. 

In the unitary operator U of the cumulant method we have to include the effect of the perturbation H\ which can 
either create or destroy pairs of spin waves. An appropriate form of U is 

U = Y[Uq, L/q = e ^« b i q -"<« b -<«). (36) 
q 

Note that the operators J7 q commute with each other because a q and 6 q are independent bosons (linear spin-wave ap- 
proximation). The coefficients /i q have to be determined by solving the equation = (<p n \ A H U\4> n ) c for appropriate 
operators A and eigenstates \<j> n ) of Hq. After transforming the cumulants into normal expectation values according 
to Appendix B we have to apply U to the eigenstates |</ , {n aq ,n i)q }) of Hq. In Appendix D it is shown that f/ q can be 
factorized into the following form: 

[/ = e (tarih/i q )a q "6i qe -(hicosh/x q )(a q "a q + 6± q 6_ q + l) e -(tanh/x q )a q 6_q ^ 
If one applies U to the unperturbed ground state |</>o) one finds 

U\<h) = (jl^/^tj e E ^ q< -|0o). (38) 

with the substitution ^ q = tanh/i q . This transformation, is the same as used within the ground-state cumulant 
formalism to describe the Heisenberg antiferromagnet, seeE3. It can be used to determine the coefficients z^ q . From 

= (Ma q 6- q )- !Z[/|0o) c (39) 
we obtain a quadratic equation for each ^ q . It has the solution0 

tanh Mq = „ q = _L(-(l + e ) + ^{l + ef - 7 (q)2) . (40) 

With these values for the coefficients the unitary transformation U (|36| ) equals the usual Bogoliubov transformation 
for the antiferromagnet which diagonalizes the Hamiltonian (|34|). Therefore we expect to obtain exactly the results 
known from linear spin-wave theory. 

In order to calculate the energy spectrum of the system we have to evaluate 

This is done in Appendix E. It leads to the following result: 

E {n a q,n b «} = E Neel + 2J ^("-aq + «6q) + X! (™«q + n 6q + X ) 

q q 

= E + ^{n a q + ™fcq) W q (42) 
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where w q is the spin-wave energy given by 



^ = y(He + 7(qK) = ~f V(l + - 7(q) 2 • (43) 

i?o denotes the ground-state energy of the Heisenberg antiferromagnet within linear spin-wave approximation. Eq. 
( [42| ) exactly equals the result of linear spin-wave theory. It describes the spectrum of a system of independent bosons 
with the dispersion w q . 



B. Weakly doped case 

Now we turn to the discussion of the weakly doped system. The t-J Hamiltonian with anisotropic magnetic exchange 
can be decomposed as follows: 

Hq ^Ising ~i~ H-Zeeman 

= J (1 + e) J2 ( S i 5 i - ^ ) + (- £ S? + £ S/), (44) 

Hi = H t + H ± 

= -* E ( & ^ + & ) + 1 E ( + ) ■ 

<ij>,a <ij> 

We have added a staggered field i?^ in ^-direction which will later be used to determine the magnetization of the 
system. 

The ground state \4>o) of H is a Neel state with M — 5-N holes where 5 is the hole concentration and N the number 
of lattice sites. The holes have fixed momenta k m and are located on the sublattice a m =T, J.. Excited eigenstates of 
H can be defined by (compare ([35]) ): 

T-r (a+) n «« (b+) nb o 

\<t>{VLa},{n a ^n b ^) = [[ , F , ■ [ [ Ck m a m \<t>N6el) (45) 

q v n aq- v nfc q- m =l 

The set of hole momenta (and spins) in this state is denoted by {kcr}. Their distribution can be different from the 
one of the ground state. The product Yl m runs over all holes in the system (k m cr m € {kcr}), and {ri aq ,?if, q } is the 
set of numbers of spin-wave excitations with momenta q on the t-and |-sublattices. 

For a proper description of the hole states we use path (or string) statesEZrt3 which lead to the concept of spin-bag 
quasiparticles. The path states are generated by repeated application of the hopping Hamiltonian H t to a hole in 
a Neel background. The moving hole leaves behind a trace of overturned spins. Here we define path concatenation 
operators A n = A n f +A„j. The operators A n f and A n ^ refer to the two sublattices. A n ^ operating on the Neel state 
with one hole, Ci^tjiNeei) , moves the hole n steps away and creates a path or string of n spin defects attached to the 
transferred hole. Explicitly, the operators A„f are defined by 



An = -7= E^i 



1 

V Zq (z 



1st = , = =, E^^/c+^%, (ier, j el, I el, m el) (46) 



^3T - i = E °™l S l S t^il R ml R lj R ji ' 

V^O (2d - 1) lj7m 

The operators A n ± for the 'down' sublattice are defined analogously with all spins reversed. Zq=4 denotes the number 

of nearest neighbor sites in the lattice. The matrices Rji and By allow the hole to jump to its four nearest neighbors 
in the first step and to only three new nearest neighbors by hopping forward in each further step: 

f 1 i, j nearest neighbors 
d% ~~ I otherwise ' ' 

( 1 j, I nearest neighbors and I ^ i 
l i I otherwise 
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The staggered magnetization per spin and /j,b (pb - Bohr magneton) has to be obtained from the free energy F by 

OF 



M, 



eff 



N ^ 1 ^ Nfi B dB A 



(48) 

B A =0 



where p denotes the statistical operator of the full system described by H. 

In the doped system we have two possible contributions to the reduction of the order parameter: the-spin bag, which 
forms the quasiparticle around each hole, and spin-wave excitations. From ground state calculationslilrEj, however, it 
is known that the main contribution to the decrease of the magnetization arises from spin waves which admix into 
the ground state. We expect that also the decrease of the magnetization with increasing temperature comes primarily 
from excited spin waves. These spin waves will be treated as bosonic excitations, i.e., spin- wave interactions will be 
neglected. Within our model we shall obtain a spin-wave spectrum which depends on the number and momentum 
distribution of the hole quasiparticles. 

Since the direct contribution from each spin-bag quasiparticle to the decrease of the magnetization is small, we 
assume it to be independent of temperature and hole concentration. This means that structure and size of the 
quasiparticle (for fixed momentum) do not change with temperature and doping. Consistent with this we employ a 
rigid-band approximation for the quasiparticles, i.e., we keep the hole dispersion ek fixed (at its value at 6 = and T = 
0). In the ground state all holes are found near the minima of the dispersion being located at k = (±7r/2, ±tt/2) (Fermi 
sea). Considering excited hole states we take into account particle-hole excitations within the lowest quasiparticle 
band. However, excitations to higher bands (interband excitations) are neglected because such excitations would have 
energies of order 2J and the temperatures we are interested in are smaller than J/2. 

In the unitary operator U of the cumulant method we have to include the effect of the perturbation Hi . Its part 
H± creates or destroys pairs of spin waves. H t provides hole hopping which can be described via the path operators 
(Eq). Within a factorization approximation U takes the form 



^holes) _ e £„ A " A ». ( 49 ) 

This ansatz generalizes the operator U from ( |36| ) to the case of weak doping. In the hole part l/( holes ) we only 
include path creation operators because we consider intraband but not interband hole excitations. We will see that 
the (intraband) hole excitations affect the spin- wave spectrum, i.e., the spin- wave renormalization depends on the 
distribution of the hole momenta in the system. Within our factorization approximation we therefore have to use 
coefficients /z q in JJ^ spm ^ which depend on the set of hole momenta {kcr}. Path operators A n are included up to length 
n = 2. The path coefficients A„ should depend on the hole momenta in ( ffBl) and on the spin- wave excitations. As 
discussed above, we expect this dependence to be weak. It is neglected wiihin our rigid-band approximation which 
fixes the values of the A n to those known from the ground-state calculation^. 

With the ansatz ( |49| ) we can now calculate the coefficients /z q for all distributions of hole momenta. In analogy to 
the undoped case we use the equation 

= (0{k<r},{O,O}|(aq&-q)' •H"t / l0{k<r},{O,O}) C (50) 

which is provided by the cumulaiiL method. The expectation values are calculated to first order in 5 using linear 
spin-wave theory, for details seeEJ'EJ. 

Having found the values for the coefficients /x q (the A„ are fixed by the rigid-band approximation for the hole 
quasiparticles) we can calculate the energies of the system states. These are described by the set of hole momenta 
{ku} and the numbers of spin- wave excitations {n aq ,rib q }. The energies are calculated in analogy to Appendix D, 
they can be written as 

E {ka},{n a ^,n b ^} = (0{k ( T},{n oq ,n i)(1 }I^C / l0{k CT },{« a<1 ,n 6 < 1 }) C 

= E NieX {1-5)- (2.7(1 -S) + B A )N + 

+ E 6k "» + ^(naq + n b<i + 1) U)q . (51) 
m q 

The renormalized spin- wave energies cD q depend on tha-JMmiber of holes M — 5 ■ N, their momenta and spins {kcr}, 
and the external field Ba- Explicitly they are given byt30 







^ = J -f{\ 5) (1 + e + 7 (q)i/ q (l - 5) ) + \\ l tF 1 {\ 8) + A 2 ~ F 2 (l - 5) 2 



2./CI ■.>)\/A3-7(q) 2 -Ai|(Af 1 -7(q)i?3) (52) 



with the substitutions tanh/i q = f q , 



A = 1 + £+ 27(T^) ^ 



A/ 



and 

iF x = Y, ( m IK L q)' ^i«6i q )'|m) 
J 



m— 1 

•4- II 



2 F 2 = ^ (m|(o q &_ q )-Jf x 4i(<*'- 

m— 1 

M 

tF 3 = J2 (m|(a q 6_ q )-ff t A 1 |m) c . (54) 



m— 1 



In the second step of ( p2| ) we have inserted the solution for 2/ q obtained from eq. (p0|). Due to the approximation of 
independent holes the expectation values with \</)q) in (|l],[32]) factorize into contributions from each of the holes. So 
the cumulant expectation values in ( [54] ) have to be taken with one- hole states \m) defined by 

I'm) = C km(7m \(f> N eel)- (55) 



The sum over m runs over all holes in the system. The expression ( |52j ) includes all terms up to first order in the hole 
concentration 5, i.e., hole- hole interactions are neglected. 

To calculate the magnetization of the system according to ( ^8|) we have to sum up the contributions from all states. 
The partition function is given by 

Z = Y Y e' l3E{ ^ i ^'^ } . (56) 

{ko-} {n a q,n bq } 

In the partition function we can exactly sum the contributions from the excited spin waves which are bosons. Note 
that the spin-wave excitations do not affect the hole part because of our rigid-band approximation. Thus we can 
define a probability p^ko-} f° r one-hole momentum distribution {krr}: 

{n a q,n bq } 

exp(-/?2: m ek m )n q (sinh^ q /2)- 2 
E {kCT y exp(-/3 £ m , e km ,) n q (sinh/to q /2)- 2 

where we have used 



(57) 



y -Kn+i)a = -J—^. (58) 
^ sinh/3w/2 v ; 

n=0 ' 1 

Finally, one obtains the following expression for the staggered magnetization (fl8|): 

Meff = \ - SM holes - £ P{ka} lW||5L cosh/3d} q /2 - l) (59) 

{kcr} q ^ A ' 

where 5M^ ; es is the contribution of the hole quasiparticles to the decrease of the magnetization. Within the rigid- 
band approximation it is independent of the temperature. It is determined by the size of the quasiparticles (spin 
bags) and can be expressed by the hole concentration 5 and the path coefficients A n of our ansatz (|49"|): 
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SM holes = S ■ ^ " . (60) 

The last term in (^9|) contributing to the magnetization is calculated numerically. The integral over the spin-wave 
momenta has to be treated carefully for small q. For an isotropic system (e = 0) and finite temperature it diverges 
logarithmically indicating an instability of the magnetic order in accordance with Mermin- Wagner theorernEi The 
sum over the hole configurations {k} has been calculated for a finite lattice. The contributions from states with hole 
excitations, i.e., states with a hole momentum distribution different from the ground state, are rather small for the 
low temperatures considered here. 

In Fig. 1 we show the staggered magnetization calculated from ([39|) for t/J — b and an anisotropy e = 2 * 10~ 6 in 
dependence on temperature and hole concentration. At T — and 5 — it has a value of about 61% of its saturation 
value which is of course identical with the result obtained in linear spin-wave theory. The small anisotropy does not 
have an essential influence on this value. For increasing temperature and/or increasing doping level the magnetization 
decreases and reaches zero at a line in the T — (5-plane. This line shown in Fig. 3 can be interpreted as the boundary 
of the antifcrromagnetic phase where we expect a second-order phase transition to a paramagnetic phase. Note 
that linear spin-wave approximation becomes questionable for vanishing sublatticc magnetization. Therefore, the 
results presented here are reasonable in the region of small doping and low temperature, but the actual calculation 
can not provide a reliable description of the system in the vicinity of the magnetic phase transition (T rs T/v) and 
beyond it. Our starting point was the Neel state, and spin fluctuations are included by a jaerturbational method 
based on cumulants, but we always describe a system with antiferromagnetic lons-range orderlij. The behavior of the 
zero-temperature magnetization is shown in Fig. 2 for e = 0, for a discussion seetll. 



V. CONCLUSION 



The aim of this work was to present a method for calculating thermodynamic properties in correlated electronic 
systems. This method is based on the introduction of cumulants and is an extension of a cumulant approach for 
ground-state calculations which has been applied in the past to a wide range of strongly correlated systems. 

After a short review of the ground-state version we have developed the extension of the cumulant method for the 
calculation of excitation energies. It is based on the introduction of an unitary operator U transforming the set 
of unperturbed eigenstates of H$ into the full eigenstates of H . For the unitary operator we use an appropriate 
exponential ansatz. Having calculated the spectrum of excitation energies we obtain the free energy and other 
thcrmodynamical quantities. The method is especially appropriate for low-temperature expansions if one has to take 
into account only low-lying excitations of the considered system. 

In the second part of this paper we have applied our method to the weakly doped t-J model with an anisotropic 
exchange at_iinj±e temperature. Based on the investigation of the renormalization of the spin-wave spectrum by 
mobile holeaiU'EU we have calculated the staggered magnetization as function of doping, temperature and magnetic 
anisotropy. These calculations are only valid within the antiferromagnetic phase of the t-J model at small doping 
concentrations and low temperatures because of the use of linear spin-wave approximation. Possible improvements 
could include magnon-magnon interaction, i.e. nonlinear spin-wave theory, which leads to a temperature-dependent 
spin-wave spectrum. 

Summarizing, we have developed a general formalism for the investigation of thermodynamic properties in both 
weakly and strongly correlated many-body systems. 
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APPENDIX A: GENERALIZED CUMULANTS 

It was proposed by KladkcEl to define cumulant expectation values with different bra and ket vectors. They are 
defined in analogy to (0) by 



ii 



The vectors ($| and \<p) must have a non-zero overlap, ($|</>) 7^ 0. It is easy to show that the cumulant equations 
( |lT| . |l3|) and ( |l^ , ^o|) are also valid with a bra vector (<i>| different from \4> n ), e.g. 

E n = ($\HU\<t> n ) c (A2) 
O = ($|(5 M -5t)t J ffC/|0 n ) c (A3) 

with (&\4> n ) 7^ 0. These equations follow from the fact that U transforms the unperturbed state \(j> n ) into an eigenstate 
of the full Hamiltonian H . 



APPENDIX B: EVALUATION OF CUMULANT EXPECTATION VALUES 



In this appendix we show how to evaluate cumulants containing an exponential ansatz for the wave operator f2 or 
the unitary operator U . The basic relation written with generalized cumulants is 



(Bl) 



with Ai and S being arbitrary operators and ( < I > |</>) ^ 0. Note that on the l.h.s. the operators S and Ai are subject 
to cumulant ordering whereas on the r.h.s. only the A, operators are cumulant entities. However, the cumulants on 
the r.h.s. are formed with the new ket vector \e S (b). 



Eq. (EUi can be proved either by integrating infinitesimal transformations (1 + S/N) and using properties of 
cumulanta3 or by explicit use of the definition of cumulant expectation values. Here we demonstrate the second way. 
Starting from the definition of generalized cumulant expectation values ( |Al[ ) for a product of arbitrary operators Ai 
and arbitrary states ($|, \cj>) with ($10) ^ we consider the following expression: 



i 



n=0 



E 

71 = 



d_ 



n 



ln($|[]e A ' A 'e« s | 



(B2) 



£=0 



The last expression can be interpreted as a series expansion with respect to £ around of the term in the brackets 



i \ i ^ 1 ' / i 

= ($| l[A? \e aS <t>r . 



|A != 0Vi 



(B3) 



In the last equation we have reintroduced generali zed cumulants, now formed with the bra state ($| and the ket state 
\e a <fi). With a = 1 we ob tain the desired result (Bl). 
Explicitly we find from (Bl): 



($| Ae 



S I A\c 



($1 AB e 



(<f>\At 



S I JAC 



($1^10) 

($| AB e s \q 



($| Ae s \<f>) ($| B e s \tf>) 
($|e s |</>) 2 



(B4) 



and so on. 
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APPENDIX C: COMPARISON OF THE ZERO-TEMPERATURE CUMULANT APPROACH WITH 

OTHER METHODS 



For practical calculations the zero-temperature cumulant method together with the exponential ansatz for the wave 
operator fl consists of selecting an appropriate set of operators S u , i.e., writing down an ansatz for the ground-state 
wavefunction. Then the coefficients a v are determined using the equations (Jl3|). The main advantage of this procedure 
compared to other methods is that the exponential term occurs only once in all equations. 

In a standard variational calculation one uses an ansatz for the wavefunction and minimizes the ground-state energy 
by variation of the coefficients. In such a calculation the ansatz wavefunction (including the exponential operator) 
usually occurs four times, Eq — (ipo\H\ipo) / (V'olV'o)- Furthermore, a wavefunction with an exponential ansatz usually 
cannot be normalized. So both numerator and denominator of the energy expression might diverge with an exponential 
of the system size whereas their ratio should be proportional to the system size. 

The physical difference between both methods is the following: In a variational calculation the aim is minimizing 
the total energy of the system whereas in the cumulant method the aim is finding an eigenstate of H . (Note that eq. 
( |l3| ) is exactly the condition of e s \4> ) being an eigenstate of H.) 

There is a close relationship of the equations ( |ill , [i2||l5| ) to the so-called coupled cluster method. This approach 
which was originally imaaited for studies in nuclear physics is also size consistent and does not involve Wick's theorem. 
For a review see BishopEd. Recently it was showntl that the coupled-cluster method can be derived from the cumulant 
expressions (|ll| ) and (^3|). Comparing practical calculations the cumulant method with an exponential ansatz is again 
easier to handle than the coupled-cluster scheme because the exponential term occurs only once in the cumulant 
equations and twice in the coupled-cluster equations. 

Usually these different methods lead to different (approximate) results when calculating ground-state quantities. 
However, if the ansatz for the ground-state wavefunction covers the exact ground state, i.e., if the subspace spanned 
by the operators S v contains the exact ground-state wavefunction, then of course all methods lead to the same exact 
result. In the following we briefly show how to derive coupled-cluster and variational equations from the cumulant 
method if one assumes that the exact ground state has the form \tp ) = e s |0 o ) with S = ^ v a v S v . We note that 
eq. (|Ts| ) also holds for arbitrary composite operators, e.g., = {4>o\ABH£l\(j)o) c for arbitrary operators A and B. 
Inserting (O) into ( |TT| ) one obtains 

E Q = {<j> \He s \<j> ) c = (0 o |e- 5 tfe 5 |0o) c = <0o|e 5 W|</> ) c ■ (CI) 
Evaluating the cumulants in analogy to appendix B leads to 

E = me- S He^ ) = ^j^. (C2) 

These are the energy expressions for the coupled-cluster and the variational scheme, respectively. The equations for 
the coefficients are obtained from (|l3|) as follows: 

= (0ol^#e s |0 o ) c = (0 o |S'+e- s J ffe s |0o) c = SlHe s \^) c . (C3) 

Transforming again the cumulants and using {<j>o\S v \^o) = one finds 

= {^o\Sle- s He s \<p ) (C4) 

and 

= {MStH\^) c + (^ \HS v \ij y = TT-^o- (C5) 
These two conditions are the equations for the coefficients a v within the coupled-cluster and the variational method. 



The second step of (C5) includes evaluating the new cumulants with |-0o) yielding exactly the four terms arising from 
the differentiation of the energy expression (^>q| HUpp) / {ipo\tpo) with respect to a„. 

We want to note here that the wave operator (n0|) of the cumulant approach is not limited to an exponential form 
(as is the case e.g. in coupled-cluster calculations). So the cumulant method appears to be the more general and 
powerful scheme for the calculation of ground-state properties. For modified applications of the cumulant approach 
see e.g.Q. 
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APPENDIX D: UNITARY TRANSFORMATION FOR THE HEISENBERG ANTIFERROMAGNET 



In this appendix we consider the Bogoliubov transformation which diagonalizcs the Hcisenbcrg antifcrromagnet 
within linear spin- wave theory. We want to prove the equivalence of the forms of Uq given in (|3^) and (37). We have 
to prove 



exp(/i(a + 6 + — ab)) 

= exp((tanh /i)a + b + ) exp (—(In cosh /i)(a + a + b + b + 1)) exp(— (tanh /x)ab) 



(Dl) 



where we used the short-hand nota tion s a and b for the boson operators a q and 6 q . 

Note that the left-hand side of (Dl) is the unitary transformation C/ q (for fixed momentum q) transforming the 
"original" bosons a and b into the "new" bosons a and (3, 



C/ q a — a — a cosh /i — b + sinh fx , 
UqbU^ — j3 — b cosh /i — a + sinh ll . 

These relations can be easily checked by expanding both sides with respect to [i. With this we obtain 
£/q ab = cosh 2 /i (ab — tanh li (a + a + b + b + 1) + (tanh 2 ll) a + b + ) J7 q . 

Another commutation relation needed is 

exp((tanh /i)a + b + ) a + a = (a + a — (tanh /i)a + 6 + ) cxp((tanh fi)a + b + ) . 

which can be derived from a well-known identity of KuboEl for arbitrary operators A and H: 



du e- {x - u)H [A, H]e~ uH 



(D2) 
(D3) 
(D4) 

(D5) 



Now we can prove (Dl). We first multiply (Dl) by exp(— (tanh fj,)ab) from the right in order to define the two functions 

(D6) 



= exp(^(a + 6 + — ab)) exp (— (tanh fi)ab) 
/ r (/z) = exp((tanh^)a + 6 + ) exp (—(In cosh /i)(a + a + b + b + 1)) . 



The next step is to compare the equations of motion of /; and f r with respect to \i. One immediately finds 

^-.h = {a+b+ -ab).h + 

CM cosh fi 

d a + b + 

—f r = — f r - f r (a + a + 6+6 + 1) tanh/i. 
OfJ, cosh [i 



(D7) 



Now we insert (D3) and (D4) into (D7) leading to: 

(a + b + - (tanh/x) (a+a + 6+6+1) + (tanh 2 (i) a + b+) fi, 



d_ 



fr 



( a +6^ 



V cosh 2 , 



(tanh fi) (a + a + b + b + 1) + (tanh 2 fi) 2a + b + j f r . 



(D8) 



Using cosh 2 fx + tanh 2 /i = 1 we see that both expressions are indeed identical. Having shown that fi(f J-) a nd f r (n) 
obey the same equation of motion and the same initial condition /;(0) = / r (0) = 1 we have proved Eq. (Dl). 



APPENDIX E: CALCULATION OF EXPECTATION VALUES FOR THE HEISENBERG 

ANTIFERROMAGNET 



Here we demonstrate the evaluation of the energy expression (41). For the denominator we find by use of ( 
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i(n a q,ni, q ) 



(^{n a ^n bcl }\U\<j) { n a ^n^}) = JJ X! V 1 



— n aq +ra 6<1 + l-2i / n 6t 

,2 i' ..2M ' ' 



(El) 



For the numerator of (Eh) one obtains 



(^{llaq^Sq} n ak ~\~ 'T-bk) I (0{ri aq ,n i , q } 1^ l^{n aq ,rn,}) 



(E2) 



and 



/ mm(n ak ,ri(, k ) 

E^(k) e 



-n ak +n bk + l — 2i 



ui-H + ^ i Kk+1)(n6k+1) ) m n . . 



mm(n aq ,nf,q) 

II E 1 - 

q^k i=0 



i+1 



n aa \ I n 



(-^q) 1 I'll ' 



After straightforward algebra we find from (E3) 



{^{n^n^Wi U\<f> {natiinbtl} ) = -^7(k)^ k (n ok + n bk + 1) (</>{„ a q,„ bq }|C/|0{„ a q,„ bq }> 



(E3) 



(E4) 



Collecting all terms leads directly to (0). 
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FIG. 1. Staggered magnetization as function of hole concentration 5 and temperature T for different values of the anisotropy 
e. t/J = 5,e = 2* 10~ 6 . 




1 2 3 4 5 

Hole Concentration (%) 



FIG. 2. Staggered magnetization as function of hole concentration 8 for different values of t/J at zero temperature. The 
magnetic anisotropy e is set to zero. 
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FIG. 3. Boundary of antiferromagnetic phase for t/J = 5 and different values of the magnetic anisotropy e. This boundary 
is given by the condition M e // = 0. 
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